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ABSTRACT 

We present observations of 25 transitions of 17 isotopologues of 9 molecules toward B335. With a 
goal of constraining chemical models of collapsing clouds, we compare our observations, along with data 
from the literature, to models of chemical abundances. The observed lines are simulated with a Monte 
Carlo code, which uses various physical models of density and velocity as a function of radius. The dust 
temperature as a function of radius is calculated self-consistently by a radiative transfer code. The gas 
temperature is then calculated at each radius, including gas-dust collisions, cosmic rays, photoelectric 
heating, and molecular cooling. The results provide the input to the Monte Carlo code. We consider 
both ad hoc step function models for chemical abundances and abundances taken from a self-consistent 
modeling of the evolution of a star-forming core. The step function models can match the observed lines 
reasonably well, but they require very unlikely combinations of radial variations in chemical abundances. 
Among the self-consistent chemical models, the observed lines are matched best by models with somewhat 
enhanced cosmic-ray ionization rates and sulfur abundances. We discuss briefly the steps needed to close 
the loop on the modeling of dust and gas, including off-center spectra of molecular lines. 

Subject headings: ISM: abundances — ISM: molecules — ISM: individual (B335) — astrochemistry 


1. INTRODUCTION 

The Bok globule, B335, is a rather round dark globule 
at a distance of about 250 pc (Tomita et al. 1979). It 
is perhaps the best case for being a collapsing protostar. 
Observations of CS and H 2 CO lines (Zhou et al. 1993; 
Choi et al. 1995) were reproduced very well with mod¬ 
els of inside-out collapse (Shu 1977). To the extent that 
such models may describe the actual density and velocity 
fields in B335, this source provides an excellent test bed 
for astrochemical models. The only remaining variables in 
modeling the lines would be the chemical abundances of 
the species in question. It is even possible to trace varia¬ 
tions in the abundance as a function of radius because the 
different parts of the line profile arise in different locations 
along the line of sight. Adding the information from the 
excitation requirements of different lines provides a probe 
of the abundance through the static envelope and into the 
collapsing core of the protostar. 

On the other hand, the depletion of molecules that is 
quite apparent in pre-protostellar cores (e.g., Caselli et al. 
2002, Lee et al. 2003) warns us that molecular lines alone 
may be misleading. In the case of B335, Shirley et al. 
(2002) found that the Shu infall model that fit the molec¬ 


ular lines (Choi et al. 1995) did not reproduce the dust 
emission. They found instead that a power law density 
model with higher densities at all radii than the best fit 
Shu model was needed to fit the dust emission. We will 
consider models more similar to the best fitting power law 
as well. 

In general, the molecular lines and dust emission have 
complementary advantages and disadvantages. The lines 
can be strongly affected by depletion that varies with ra¬ 
dius, while the dust shows no convincing evidence so far 
for variation of opacities with radius (Shirley et al. 2002). 
On the other hand, variation in opacities with radius is 
also not ruled out, and the actual value of the opacity at 
long wavelengths is quite uncertain, by factors of at least 
3 and possibly more. The dust emission is sensitive only 
to the column density along a line of sight, while the line 
emission can in principle probe the volume density via ex¬ 
citation analysis. Finally, only the lines can probe the 
kinematics, but that probe can be confused by depletion 
effects (Rawlings & Yates 2001), and the dust is needed 
to constrain these effects. Clearly, the best approach is a 
unified model for both gas and dust components. 

We will present new observations of a large number of 
species toward B335, using Haystack Observatory and the 
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Caltech Submillimeter Observatory. We will also present 
the results of detailed models of radiative transport in 
dust to determine the dust temperature for several dif¬ 
ferent physical models. Next, we will calculate the gas ki¬ 
netic temperature, including gas-dust interactions, cosmic 
rays, and photoelectric heating. With these as a basis, 
we will calculate the molecular excitation and radiative 
transport, using a Monte Carlo code (Choi et al. 1995). A 
telescope simulation code will produce model line profiles, 
given an input model of the density, temperature, veloc¬ 
ity, and abundances as a function of radius, for comparison 
with the observed line profiles. Based on the comparison, 
the abundances of various species will be constrained. We 
will use step function models for the abundances and also 
the results of new calculations of abundances in a cloud 
collapsing according to the Shu picture (Lee et al. 2004). 

2. OBSERVATIONS 

We obtained observations of the HCO + and N 2 H + 
J = 1 —> 0 lines at the Haystack Observatory in 1995 
March. Observations of a large number of lines were ob¬ 
tained at the Caltech Submillimeter Observatory in the 
period 1995 March to 2001 July. Table 1 provides the ref¬ 
erence frequency for the line, the telescope, the main beam 
efficiency (r? m b)) the full width at half maximum beam size 
(6b), the velocity resolution (Sv), and the date of observa¬ 
tion. We also provide this information for several obser¬ 
vations obtained previously that are used to constrain the 
modeling. The frequencies in Table 1 are either those used 
during observing or those used later to shift the observed 
data to an improved rest frequency. For most lines with 
hyperfine components, these are the reference frequencies 
suitable for a list of hyperfine components that were used 
to fit lines. In the case of the N 2 H + J = 1 —> 0 line, it 
is the frequency of the isolated hyperfine component, best 
suited for determining the velocity. 

In the following sections, we assume that the centroid 
of B335 is at a = 19 ft 34 m 35.4 s ; <5 = 07°27'24" in 1950 
coordinates. This position agrees within 1" with the cen¬ 
troid of the submillimeter emission mapped with SCUBA 
(Shirley et al. 2000). This position was originally based 
on the position of the millimeter continuum source seen 
by Chandler & Sargent (1993); more recent interferomet¬ 
ric data find a compact component located 3"6 west and 
1" south of this position (Wilner et al. 2000). At this posi¬ 
tion, a continuum source is also seen at 3.6 cm, attributed 
to a time variable radio jet elongated along the outflow 
axis (Reipurth et al. 2002). The difference between our 
position and the position of the compact component is not 
significant for the resolution of these observations. Some of 
our data were obtained before we settled on this position. 
In cases where we have a map, we may have resampled the 
data spatially to synthesize a spectrum at the submillime¬ 
ter centroid position, resulting in a slight degradation of 
the spatial resolution. 

3. RESULTS 

The primary observational results are presented in Ta¬ 
ble 2 and Figure 1 to Figure 7. The table gives the inte¬ 
grated intensity (fT A dv), the peak antenna temperature 
(TjJ), the velocity with respect to the local standard of 
rest ( vlsr ), and the linewidth (FWHM), Av. For simple, 


single-peaked lines, these were determined from a Gaus¬ 
sian fit. For self-reversed lines without hyperfine structure 
(HCO+ J = 1 —> 0, J = 3 —> 2, and J =4 —> 3), fT\dv is 
the total area under the full line, T\ is the strength of the 
stronger peak, Vlsr is the velocity of the dip, determined 
by eye, and Av is J T A dv divided by T\. For lines with 
hyperfine structure (C 17 0 J = 2 —» 1, N 2 D + J = 3 —» 2), 
f T A dv gives the area under all the hyperfine components, 
T\ gives the peak of the strongest, usually blended com¬ 
ponents, and Vlsr and Ac come from a fit with all the 
hyperfine components. For the most complex situation, 
lines that are self-reversed, with hyperfine structure, vari¬ 
ous strategies were adopted. For CN J = 2 —► 1, / T A dv, 
T' a . and Vlsr were determined as for double peaked lines, 
but Av was determined from an isolated component. The 
spectrum of CN (Figure 6) clearly shows that the main hy¬ 
perfine line is self-reversed. For N 2 H + J = 1 —> 0, all line 
parameters were determined from the isolated component 
at the frequency given in Table 1, as suggested by Lee, 
Myers, & Tafalla (2001). 

The observations that are compared to full models are 
shown as solid lines in Figures 1 to 6. The CN J = 2 —> 1 
spectrum in Figure 6 has not been modeled, and the 
dashed line is just a fit to the hyperfine components. Other 
spectra that are not modeled in detail are shown in Figure 
7. These include spectra with complex hyperfine splitting 
that we cannot model in detail and molecules without good 
collision rates. 

The HCN J = 3 —> 2 line is peculiar in that there is es¬ 
sentially no emission at velocities that would normally be 
associated with the red part of the main hyperfine compo¬ 
nent. To ensure that this effect was not caused by emission 
in the off position (10' west), we took a deep integration 
in the off position. No emission was seen at a level of 0.08 
K. 

Single-peaked lines without overlapping hyperfine com¬ 
ponents provide the best measure of the rest velocity of 
the cloud. Based on those lines least likely to be optically 
thick, the cloud velocity is ( Vlsr) = 8.30 ± 0.05 km s -1 . 
For self-reversed lines, the mean velocity of the dip (de¬ 
termined by eye) is (vdi P ) = 8.41 ± 0.06 km s -1 . All the 
values for v c n p exceed those for ( Vlsr ), by amounts rang¬ 
ing from 0.05 km s -1 to 0.25 km s^ 1 . The mean shift, 
(vdip — (Vlsr )) = 0.11 ± 0.07. Also, lines from higher 
J levels have higher Vdip than those from lower J levels, 
suggesting that the dip arises partially from inflowing gas. 
The three lines of HCO + , for example, have their dip at in¬ 
creasing velocity, with the J = 4 —> 3 showing Vdip = 8.55 
km s -1 . This progression is similar to a pattern seen in 
CS lines toward IRAM04191 by Belloche et al. (2002). 

4. THE MODELING PROCEDURE 

We use the extensive observations described above to 
test models of the source. All the models are spherical 
models with smooth (non-clumpy) density distributions. 
We focus on inside-out collapse models, though we discuss 
some variations on this basic model. All models include 
self-consistent calculations of the dust and gas tempera¬ 
ture distributions (§4.1) and calculations of the molecular 
populations, radiative transport, and line formation (§4.3). 
Two kinds of models of the abundances as a function of ra¬ 
dius are used: step function models, and abundances from 
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an evolutionary chemical calculation (Lee et al. 2004), as 
described in §4.2. 

4.1. Determining Temperatures 

The first step in comparing a physical model to observa¬ 
tions is to determine the temperatures that correspond to 
a particular density distribution. The dust temperatures 
can be calculated self-consistently for a particular density 
distribution by various radiative transfer codes. We used 
the code of Egan et al. (1988) and the techniques described 
by Shirley et al. (2002) for constraining parameters. 

We assumed that dust opacities are given by column 
5 of the table in Ossenkopf & Henning (1994), known as 
OH5 opacities, because these have been shown to match 
many observations of star forming cores (e.g., Shirley et al. 
2002). One difference between the models by Shirley et al. 
and the current work is in the treatment of the interstellar 
radiation field (ISRF). In the previous work, we decreased 
the strength of the ISRF by a constant factor ( sjsrf ) at 
all wavelengths (except for the contribution of the cosmic 
microwave background). For B335, we used sisrf = 0.3. 
In the present work, we instead attenuate the ISRF us¬ 
ing the Draine & Lee (1984) extinction law and assuming 
Ay = 1.3 mag. This procedure affects short wavelengths 
much more than long wavelengths, leading to a somewhat 
less pronounced rise in dust temperature toward the out¬ 
side of the cloud. The choice of Ay = 1.3 mag is somewhat 
arbitrary, but it accounts for the fact that molecules re¬ 
quire some dust shielding. It will also produce consistent 
results when we consider the gas energetics. 

Shirley et al. assumed an outer radius of 60,000 AU for 
most B335 models. We will mostly use an outer radius of 
0.15 pc (31,000 AU), as used by Choi et al. (1995). Studies 
of the extinction as a function of impact parameter from 
HST/NICMOS data are consistent with an outer radius of 
about this size (Harvey et al. 2001), in the sense that the 
extinction decrease with radius blends into the noise at 
that radius. The choice of outer radius makes little differ¬ 
ence in most models. The inner radius for the dust models 
is taken to be 10 -3 of the outer radius for the dust mod¬ 
els in order to capture the conversion of short-wavelength 
radiation from the forming star and disk to longer wave¬ 
length radiation. The stellar temperature is set to 6000 
K, but this choice is completely irrelevant because of the 
rapid conversion to longer wavelength radiation. The lu¬ 
minosity was set to 4.5 L 0 , which provided the best fit to 
a Shu model in Shirley et al. (2002). 

Once we have a dust temperature distribution, Td(r), 
and a density distribution, n(r), we can compute the gas 
temperature distribution, Tk(t). This was done with a 
gas-dust energetics code written by S. Doty [see Doty & 
Neufeld (1997) and the appendix in Young et al. (2004) for 
descriptions]. This code includes energy transfer between 
gas and dust, heating by cosmic rays and the photoelectric 
effect with PAHs, and molecular cooling. 

The gas-dust energy transfer via collisions depends on 
the total grain cross section per baryon, averaged over the 
distribution of grain sizes (E<j). Following the discussion 
in the Appendix of Young et al. (2004), we take this value 
to be 6.09 x 10~ 22 cm 2 . The cosmic ray heating depends 
on the cosmic ray ionization rate (£), which we take to be 
3 x 10 -17 s _1 (van der Tak & van Dishoeck 2000). The 


photoelectric heating follows the equation of Bakes & Tie- 
lens (1994), which includes heating from the photoelec¬ 
tric effect on very small grains. The rate depends on the 
strength of the ultraviolet portion of the ISRF and the elec¬ 
tron density. Because this heating is only important on the 
outside of the cloud, we set the electron density to 1 x 10~ 3 
cm -3 . The radiation field is assumed to be attenuated by 
the surrounding medium according to tuv = 1.8Ay; with 
Ay = 1.3, the scale factor for the ISRF impinging on our 
model’s outer radius is Go = 0.1. Once inside the cloud, 
the radiation is attenuated according to a fit to the atten¬ 
uation produced by the dust assumed to be in the cloud. 

The result of these calculations is shown for a typical 
model in Fig. 8. For small radii, Tk ~ Td , as is usually 
assumed, but Tk falls below Td with increasing radius, 
as the density becomes too low for collisions with dust 
to maintain the kinetic temperature at the dust temper¬ 
ature. Then, at some radius, Tk rises as photoelectric 
heating takes over, and Tk > Td . The downturn in Tk 
at the cloud edge appears to be real and caused by the 
cooling lines of CO becoming optically thin (see Young 
et al. 2004). However, this drop in temperature has no 
appreciable effect on the resulting line profiles. 

The amount of photoelectric heating is the least cer¬ 
tain of these inputs, as the external attenuation has a 
large effect on how warm the outer cloud gets. To con¬ 
strain Go, we modeled the lower three lines of CO and 
compared to data in the literature (e.g., Goldsmith et al. 
1984, Langer, Frerking, & Wilson 1986). To avoid pro¬ 
ducing a CO J = 1 —> 0 line that exceeded the observa¬ 
tions, Go definitely needed to be decreased from unity. 
The value of Go = 0.1 provided the best match, and this 
was actually used to constrain the external extinction to 
Ay = 1.3. Changes by a factor of 2 in Go (or ±0.4 in Ay) 
produced CO lines that differed from the observations by 
about 30%, while having no appreciable effect on the lines 
of other species. While other variables are uncertain in 
the photoelectric heating, the attenuation of the ultravi¬ 
olet radiation from the ISRF is the most important vari¬ 
able; comparison to observations of CO readily constrain 
it. The results are reasonable; one does not expect sig¬ 
nificant molecular gas for Ay < 1 (van Dishoeck & Black 
1988). One could trade off the value of sisrf and the 
external extinction, as long as the effective Go is not too 
different from 0.1. Constraining these separately is diffi¬ 
cult (Shirley et al. 2004) and not particularly relevant for 
this paper. 

The cooling rates (Doty & Neufeld 1997) depend pri¬ 
marily on the CO abundance and the width of the lines 
(through trapping); we assume X[CO) = 7.4 x 10 -5 and 
b = 0.12 km s _1 , except for some tests described below. 
We note in passing the dangers of simplistic interpretation 
of observed CO lines; turning the observations into a ki¬ 
netic temperature would lead one to conclude that Tk is 
constant within the cloud, while it clearly is not. 

The parameters that describe the standard physical 
model are summarized in Table 3. 

4.2. Chemical Modeling 

Two kinds of abundance models are employed. The first 
is strictly ad hoc, using a step function to describe the 
abundance of each species as a function of radius. These 
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models have three free parameters per species: X , Tdep, 
and fdep- The abundance in the outer parts of the cloud 
( X ) is assumed to decrease inside a depletion radius ( fdep ) 
by a factor {fdep)- 

The second are true chemical models, based on the cal¬ 
culations presented by Lee et al. (2004). These calcula¬ 
tions follow the chemical evolution through an evolution¬ 
ary sequence that includes at each step a self-consistent 
calculation of the dust and gas temperatures, using the 
techniques described in the previous section. The evolu¬ 
tionary model assumes a slow build-up in central density, 
via a sequence of Bonnor-Ebert spheres, to the point of 
a singular isothermal sphere, at which point, an inside- 
out collapse (Shu 1977) is initiated. After this point, the 
chemistry is calculated for each of 512 gas parcels as it 
falls into the central region. Thus, gas inside the infall 
radius carries some memory of the conditions from farther 
out. We adopt the model of Young and Evans (2004) for 
the evolution of luminosity in order to calculate the evo¬ 
lution of dust temperature. Physical parameters in the 
model are selected to have a total internal luminosity and 
a dust temperature profile similar to those obtained from 
the dust modeling of B335, at the time step of r ln f = 0.03 
pc. The model core is assumed to stay in the same envi¬ 
ronment through its evolution with the same Ay and Go 
calculated in the previous section. The chemical calcula¬ 
tion includes the interaction between gas and dust grains 
as well as gas-phase reactions, but the surface chemistry 
is not considered in the calculation. For details of the 
chemical evolution model, refer to Lee et al. (2004). 

For both types of models, isotope ratios were con¬ 
strained so that the abundance was only a free parame¬ 
ter for the whole of the isotope complex. DCO + was the 
exception, as it is subject to large fractionation effects. 
Assumed isotope ratios are the same as those used by 
Jprgensen et al. (2004) and are given in Table 4. Wouter- 
loot et al. (2004) have recently suggested a slightly higher 
ratio for 18 0/ 17 0 of 4.1, but this value would fit our data 
on the CO isotopologues somewhat worse. 

4.3. Modeling of Line Profiles 

The line profiles were modeled with a Monte Carlo code 
(me) to calculate the excitation of the energy levels and a 
virtual telescope program (vt) to integrate along the line of 
sight, convolve with a beam, and match the velocity reso¬ 
lution, spatial resolution, and main beam efficiency of the 
observations (Choi et al. 1995). All lines were assumed 
to be centered at 8.30 km s -1 , based on the average of 
optically thin lines. 

The input physical conditions (density, temperature, 
and velocity fields) were taken from the physical model 
being tested, using the results of the gas-dust energetics 
code for Tk[t). Models require input data about each 
molecule, as well as about the source. For CS, we used 
collision rates from Turner et al. (1992). For HCO + and 
N 2 H+, we used collision rates supplied by B. Turner, based 
on his extension of previously calculated rates to higher 
temperatures and energy levels (Turner 1995). For H 2 CO 
and para-H 2 CO, we used rates computed by Green (1991). 
Rates for HCN came from Green & Thaddeus (1974) and 
those for CO from Flower & Launay (1985). In some cases, 
rates have been extrapolated to lower temperatures. 


For C 17 0, HCN, and N 2 H + , the lines have hyperhne 
structure that is partially resolved. For these lines, me and 
vt models were run separately for each clearly resolved hy¬ 
perhne component, with abundances adjusted to simulate 
the fraction of the transition probability in that compo¬ 
nent; the results were added to make the final simulated 
line. Components separated by less than the 1/e width of 
the velocity dispersion were aggregated into a single com¬ 
ponent; the aggregated components are listed in Table 5. 
This procedure captures the essence of the hyperhne split¬ 
ting, but it is not rigorous because trapping is not handled 
correctly when there is partial line overlap (see Keto et al. 
2004). 

All models were run with 40 shells. The inner radius 
was 2 x 10 -3 pc, corresponding to T/7 at a distance of 
250 pc. This radius is larger than the inner radius for the 
dust models because the molecular lines are not sensitive 
to emission from very small scales because of beam dilu¬ 
tion. The convergence criterion for populations was set 
to 10% for finding the region of best-fitting parameters. 
Final models were run with a 2% convergence criterion 
to ensure accuracy; differences between these models and 
those run with 10% accuracy were small. The minimum 
fractional population tested for convergence was ICC 6 . For 
an explanation of these criteria, see the Appendix in Choi 
et al. (1995). 

5. INSIDE-OUT COLLAPSE MODELS 

The physical properties of the standard model are given 
in Table 3. The standard physical model is the inside- 
out collapse model (Shu 1977) that best matched (Choi 
et al. 1995) the CS and H 2 CO data taken by Zhou et al. 
(1993). Choi et al. (1995) modeled CS and H 2 CO lines 
from the IRAM telescope, assuming constant abundances, 
and found a best fit ri„/ = 0.03 pc. This was a compro¬ 
mise, as H 2 CO favored smaller Ti n f than did CS. 

While the CS data were still well matched with constant 
abundance, the new data on more H 2 CO lines suggested 
enhanced abundances on small scales, as did the HCO + 
data. As a result, we tested step function abundance mod¬ 
els. 

5.1. Step Function Abundances 

To avoid too many free parameters, we required that 
fdep = finf ■ While this particular choice has no theoretical 
justification, it leaves only two free parameters per species. 
With the constraints on isotope ratios in Table 4, we are 
left with 15 free parameters for 8 species, including the spe¬ 
cial case of DCO + , explained below. The abundances in 
Table 6 are those that fit the current data reasonably well, 
as judged by eye and statistical measures. We calculated 
both the reduced chi-squared () and the absolute devi¬ 
ation {AD = | T^{model; i) — T^{obs;i)\/N) over the 
line profiles. The absolute deviation is more influenced by 
strong lines for which the shape is important to match, so 
we use it primarily, though the \r criterion does not differ 
in the choice of best model. We have not run a complete 
grid of models; instead, we employed some judgment to lo¬ 
cate regions of parameter space with decent fits to the line 
profiles. Once reasonably good fits were obtained, both 
X and fdep were varied by factors of 3 in each direction, 
showing substantially worse fits. These parameters should 



B335: A Laboratory for Astrochemistry 


5 


be considered constrained at that level. The abundances 
in Table 6 are those for the best-fitting step functions, and 
the predicted line profiles are shown in Figures 1 to 6 as 
gray lines. The values of AD in Table 6 are an average 
over all the lines for all isotopologues of that species. 

The CS lines were still matched best with constant 
abundance. On the small scales (0.003 pc) probed by inter¬ 
ferometers, the CS is clearly depleted in the envelope and 
the J = 5 —► 4 emission arises from a clump that is offset 
from the central source (Wilner et al. 2000). A model with 
CS depleted by a factor of 10 for rd ep = 0.003 pc showed 
no appreciable effect on the J = 2 —>lorJ = 3-^2 lines, 
but it did decrease the predicted J = 5 —> 4 intensity 
slightly. 

The H 2 CO abundance that best fits the data is slightly 
higher than was found by Choi et al. (1995), mostly to 
improve the fit to the lines of H 2 13 CO and para-H 2 13 CO. 
We also increased the abundance of both H 2 CO and para- 
112 CO in the inner parts of the cloud to improve the fit to 
our new CSO observations of the higher-J lines, whereas 
Choi et al. (1995) had found a constant abundance to 
be satisfactory. Even so, we do not reproduce the very 
high excitation H 2 CO lines, indicating that a warm, dense 
region must exist that is not predicted by the basic model. 

The abundance of H 2 CO listed in Table 6 is actually 
the abundance of ortho-H^CO. Minh et al. (1995) found 
that ortho-H 2 CO/para-H 2 CO was 1.7 in B335, assuming a 
uniform cloud. Our modeling, which employs density, tem¬ 
perature, and velocity gradients, confirms that this ratio 
works well in reproducing the observations, but we have 
not determined the range of acceptable values. Minh et 
al. (1995) noted that this ratio was consistent with ortho¬ 
para equilibration on cold dust grains and suggested that 
the gas-phase H 2 CO in B335 had formerly resided on dust 
grains. In this picture, they suggested that warming by the 
newly-formed star or by shocks had liberated the H 2 CO 
from the dust grains. Our model for the dust temperature 
indicates that Td stays below 20 K until r < 0.006 pc (6"), 
where any H 2 CO would be beam diluted. Thus, other 
means for releasing the H 2 CO from dust mantles should 
be explored. 

The lines of HCO + and its isotopologues are best 
matched with a model with increased abundances in¬ 
side rinf. HCO + is clearly quite abundant in B335, as 
witnessed by the detection of HC ls O + . The observed 
J = 3 —► 2 and J = 4 —> 3 lines of HCO + are somewhat 
weaker than the models predict and the dip is shifted to 
the red (§5.2). The abundance of DCO + was treated as 
a free parameter but rdep and fdep were constrained to 
the same value as for HCO + ; the best fit was obtained for 
HCO+/DCO+ of 55. 

The C 18 0 lines were fitted best with decreased abun¬ 
dances (but only by a factor of 3) inside r ln f. With the 
enforced C 18 0/C 17 0 ratio of 3.5, the C 17 0 model line 
is a bit weaker than the observed line, but the data are 
rather noisy. Using the standard isotope ratio, the best 
fit abundance of C ls O would imply X(CO) = 4 x 10~ 5 . 
This abundance is substantially less than expected from 
chemical models and even less than what we assume in 
our calculation of cooling rates (§4.1). To see the conse¬ 
quences, we ran a model with the abundance adjusted to 
this value in the calculation of cooling rates. The value of 


Tk in the outer parts of the cloud was increased by a few 
degrees, but the effect on most molecular lines was very 
small, indicating that the best fit is not affected by this 
slight inconsistency. The CO line predictions were excep¬ 
tions, as these lines actually got stronger with decreased 
CO abundance because of the higher Tk in the relevant 
layers of the cloud. 

The two lines of N 2 H + both have hyperfine structure 
and our method for dealing with this is only approximate. 
Nonetheless, the J = 1 —» 0 transition is matched reason¬ 
ably well (Fig. 6) with a factor of 10 increase in abundance 
inside ri n f. In contrast, the satellite hyperfine lines of the 
J = 3 —> 2 line are clearly stronger than the models can 
explain. 

The most troublesome species was HCN. The satellite 
hyperfine lines of both J = 1 —> 0 and J = 3 —> 2 transi¬ 
tions are much stronger than the models can account for, 
even with a very large HCN abundance. Still larger abun¬ 
dances predicted lines of the stronger hyperfine compo¬ 
nents that were much stronger than observed. In addition, 
the HCN J = 3 —> 2 line is very peculiar, with the red side 
of the line essentially missing, indicative of a deep absorp¬ 
tion layer. To try to match some of these features within 
the constraints of our model, we depleted HCN by a factor 
of 10 inside r,„/. This helped, but the fits are still poor. 
The fact that the H 13 CN J = 3 —> 2 line was not detected 
makes the strength of the hyperfine satellite lines (Fig. 7) 
even harder to understand. 

5.2. Variations in the Physical Model 

With the additional freedom of the step function abun¬ 
dance profile, is r, n / = 0.03 pc still the best model? This 
question was explored to a limited degree; for each new 
ri n f, the abundances of each species were optimized, but 
the shape of the step function was not allowed to change, 
except for CS, where changes in fdep were allowed. For 
modest changes (rinf = 0.02 — 0.04 pc), the overall fits 
were not much worse. As found by Choi et al. (1995), the 
CS favored n n f = 0.03 pc, while H 2 CO favored smaller 
ri n f. For factor of 3 changes, the fit degraded substan¬ 
tially (Fig. 9). For r,; n / = 0.01 pc, optically thick lines 
were too narrow and the two peaks were nearly equal in 
strength, unlike the observations. For r, n / = 0.09 pc, 
those lines were too wide and the blue/red ratio was too 
large. There was also a greater conflict between the re¬ 
quirements of optically thick and optically thin lines; if the 
abundance was increased to match the latter, the former 
became too strong. Within the constraints on abundances 
that we imposed, infall radii different by a factor of 3 would 
be strongly ruled out. The mean absolute deviations over 
all species ((AD)) for these different models are listed in 
Table 7. 

The constraints on the infall radius from the molecu¬ 
lar line observations are inconsistent with those found by 
modeling the continuum emission (Shirley et al. 2002). 
The predicted intensity profiles of the model from Choi et 
al. (1995) were too flat to match the observations at 850 
and 450 /im (see Fig. 6 and Table 3 of Shirley et al 2002). 
To make an inside-out collapse model fit the data, Shirley 
et al. had to use a very small infall radius, r = 0.0048 
pc, more than 6 times smaller than the infall radius that 
matches the line profiles. Our modeling confirms that this 
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small infall radius cannot match the line profiles. This 
fundamental discrepancy between the models of the dust 
and molecular line emission will be discussed further in §6. 

Harvey et al. (2001) found that an inside-out collapse 
model with rj„/ ~ 0.03 pc fit the extinction data well, but 
only if the density was increased everywhere by a factor of 
about 5. We tried this model; decreases in abundances by 
about an order of magnitude were required for most species 
to bring line strengths back to near observed values. We 
had to decrease the CS abundance within rdep to match 
the data, but this change is not unreasonable. The aver¬ 
age deviation is somewhat larger than the standard model, 
but not terrible. Without constraints on abundances, it is 
hard to rule out variations in the physical model of this 
magnitude. However, the shapes of the CS lines were not 
reproduced well (top panel in Fig. 9), with blue/red ratios 
clearly less than the observations. 

Moving farther afield, one may consider other collapse 
models. In some sense, the opposite extreme to the Shu 
(1977) model is the Larson-Penston similarity solution 
(Larson 1969, Penston 1969). Line profiles from this model 
were generated by Zhou (1992) and found to be consider¬ 
ably wider than those observed in regions of low-mass star 
formation. More recently, Masunaga et al. (1998) have 
shown that radiation hydrodynamical (RHD) calculations 
of collapse are well approximated by a modified Larson- 
Penston model. This model produces lower infall velocities 
than the original Larson-Penston model, which decrease 
with radius. Masunaga & Inutsuka (2000) have simulated 
line profiles from the RHD models, finding blue profiles 
and smaller linewidths, qualitatively consistent with those 
seen in low mass cores. While these models may indeed 
have application in some regions, the linewidths listed in 
Table 4 of their paper for models after formation of the 
central core are larger than those in B335 by factors of at 
least two. 

There are hints in the spectra of deviations from the 
Shu model, particularly in the shift of Vdi P to higher ve¬ 
locity for lines of higher excitation. This shift can be see 
more clearly in Figure 10, where three lines of HCO + are 
shown. The best step function model is also shown; it 
does not reproduce this shift. The dip is caused by ab¬ 
sorption from low-excitation material. In the Shu model, 
the outer, static envelope is dominating this absorption. A 
model with inward motion in this outer layer might better 
reproduce this shift. 

5.3. Self-consistent Chemical Models 

The chemical models are constrained by assuming an 
entire evolutionary history for the core, as detailed by Lee 
et al. (2004). We consider only their standard model 
of the evolution of physical conditions and luminosity to 
define the physical conditions, including the dust temper¬ 
ature profile, at the time step for which = 0.03 pc. 
Compared to the standard model of Lee et al., we allowed 
adjustment of only 3 free parameters: the binding energy 
to the dust, set by the assumed nature of the dust surface; 
the initial abundance of elemental sulfur; and the cosmic- 
ray ionization rate. The different models are summarized 
in Table 8. 

First, binding energies of molecules onto three different 
dust grain surfaces were checked. For this comparison, the 


initial elemental abundances and the cosnric-ray ionization 
rate were the same as those in the standard model of Lee 
et al. (2004). For the CS lines, the binding energy onto a 
CO-dominant grain mantle works the best and the value 
of (AD) is slightly better than for Si0 2 . However, the 
low binding energy of molecules onto the CO mantle leads 
to less freeze-out of CO, and, in turn, N 2 H + is destroyed 
by abundant CO in the gas phase. As a result, simulated 
CO isotopologue lines are too strong, and simulated N 2 H + 
lines are too weak compared to the observations. Attempts 
to improve the fit by reducing the initial abundance of car¬ 
bon make the fit to CS worse while still not making the 
models fit CO and N 2 H + profiles. 

At the other extreme, CO and CS are frozen-out signifi¬ 
cantly onto H 2 O-dominant grain mantles to produce much 
weaker lines than the observations indicate. HCO + , which 
is a daughter molecule of CO, is also depleted from the 
gas phase. Although N 2 H + is less likely to be destroyed 
by CO, even nitrogen molecules are easily frozen-out on 
the H 2 0 mantle, decreasing the N 2 H + abundance. In ad¬ 
dition, HCN increases by 3 orders of magnitude at radii 
smaller than 0.004 pc, compared to the abundance in the 
outer regions, giving a very broad line wing, which is not 
present in the observed lines. Except for weaker CS lines, 
the lines simulated with abundance profiles calculated for 
bare Si0 2 grain surfaces show much better fits to actual 
data than do those from other assumptions about grain 
surfaces. 

We adopt the bare Si0 2 grain surface as our standard. 
In this model, the CO evaporation radius is about 0.006 
pc. At radii less than this radius, almost all CO is des¬ 
orbed from dust grain surfaces. Next, we varied the initial 
abundance of sulfur to improve the fit to the CS lines. An 
increase of the initial sulfur abundance by a factor of 5 
gives the best results with the Si0 2 grain surfaces. Other 
molecular lines do not vary much with the initial abun¬ 
dance of sulfur. The abundance profiles of this model are 
shown as solid lines in Figure 11. In all chemical models, 
the N 2 H + and HCO + lines are weaker than the observed 
lines, so we tested various cosnric-ray ionization rates. A 
cosnric-ray ionization rate increased by a factor of 2 pro¬ 
duced the best fit for N 2 H + and HCO + lines. We increased 
the ionization rate in the energetics calculation for consis¬ 
tency. HCN lines become stronger with the cosmic-ray 
ionization rate, and (AD) becomes somewhat worse if we 
increase the ionization rate by a factor of 5. The abun¬ 
dances for the model with ionization enhanced by a factor 
of 2 are shown by the dotted lines in Figure 11. 

The chemical models produce abundances with large 
(many orders of magnitude) variations with radius and 
quite complex radial structure. For explanations of these 
effects, see Lee et al. (2004). Note in particular the large 
decreases in abundance at large radii for most species 
caused by photo-dissociation. The large decrease in CO 
abundance over a wide range of radii reflects freeze-out 
onto grain surfaces, and some other species follow this 
trend, but N 2 H + behaves oppositely because CO destroys 
N 2 H + . Likewise, most species show a peak at small radii, 
where CO evaporates, because those species also evaporate 
there, while N 2 H + decreases when CO evaporates. The 
abundances for the step function model are also shown in 
Figure 11. In some cases, they are dramatically different. 
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The line profiles resulting from the best fit chemical 
model are shown in Figures 1 to 6 as dotted lines. The 
average values of absolute deviation for the best chemical 
model are also listed in Table 6 for comparison to those 
from the step function models. In most cases, the values of 
AD are slightly worse, most notably for DCO + and H 2 CO, 
but many are similar to those of the step function models. 
Compared to the best fit with step function abundances, 
the best fit with chemical abundances shows less deep ab¬ 
sorption dips in the CS lines, weaker lines of the higher 
transitions in N 2 H + and H 2 CO, and less blue asymmetry 
in HCO + lines. In addition, the predicted lines of HCO + 
isotopologues are weaker than in the step function models. 
To match the observations, the isotope ratios would need 
to be increased by factors of 3 to 5. The deuterium ratio 
for the DCO + 3—2 would also need to be increased by a 
factor of 5 to match the observed line. These discrepancies 
result from the fact that B335 has quite strong lines of rare 
isotopologues of HCO + . However, standard isotope ratios 
for C 18 0, 13 CO, and C 17 0, or H 2 CO and H 2 13 CO produce 
reasonably good matches to the observations. Chemical 
abundances predict still weaker H 2 CO lines in high exci¬ 
tation transitions than do step function abundances. As 
mentioned in §5.1, higher densities at small radii are nec¬ 
essary to account for the weak high excitation H 2 CO lines. 
Also, the N 2 H + 3—2 line is relatively weaker than the 1—0 
line when compared to the observed lines, again suggesting 
the presence of higher densities at small radii. HCO + and 
HCN lines simulated with chemical abundances are nar¬ 
rower than the observed lines. These two molecules are 
abundant in outflowing gas. Therefore, the lines might 
be affected by the outflow, which is not considered in this 
work. The satellite hyperfme lines of HCN 1—0 are pro¬ 
duced better with the abundance profile from the chemical 
model than with a step function. However, even chemical 
abundances cannot predict the satellite hyperfme lines of 
the HCN 3—2 as well as the absence of the red compo¬ 
nent of the main group. This result suggests the existence 
of a region with very high HCN abundance. The chem¬ 
ical models also do not reproduce the shift of Vdip with 
increasing J seen in the HCO + (Fig. 10). 

We also tested different time steps in the same lumi¬ 
nosity model. In time steps earlier than the time step 
for r vn f = 0.03 pc, the model core has higher densities 
at small radii, and the total internal luminosity is smaller 
than observed for B335 and vice versa for Ti n f > 0.03 pc. 
According to the test, H 2 CO and N 2 H + lines are fitted 
much better with the chemical abundances in an earlier 
time step for n n f = 0.015 pc. The satellite hyperfme lines 
of the HCN 3—2 are also well fitted, while its main group 
is too strong. Higher densities at small radii cause these 
results. In addition, at this time, the CO evaporation ra¬ 
dius is about 0.004 pc, so less CO evaporates compared 
to the time step for r in f = 0.03 pc. As a result, N 2 H + is 
more abundant. However, at this time, the infall veloci¬ 
ties are not big enough to produce the degree of the blue 
asymmetry in CS lines. The CS lines predicted by models 
with three times smaller and larger are also shown 
in Figure 9 and the values of (AD) are given in Table 7; 
unlike the step function models, the luminosity is different 
for these other values of rj„/, because the luminosity in¬ 
creases with time. Nonetheless, similar problems to those 


encountered in the step function models appear at earlier 
and later times. 

6. DISCUSSION 

Both step function models and evolutionary chemical 
models do a reasonable job of fitting most of the data. 
Models with constant abundances are not adequate for 
fitting most observations. In addition, the evolutionary 
chemical models clearly indicate that abundances vary by 
orders of magnitude as a result of freeze-out and desorp¬ 
tion (Lee et al. 2004). These conclusions are similar to 
those of Jprgensen et al. (2004), who find that constant 
abundance models are unsatisfactory and that a drop func¬ 
tion works better. The drop function, though simpler, is 
similar in shape to the abundance profile of CO in Figure 
11. It allows a region of lower abundance at intermediate 
radii and a return to high abundances at small radii. The 
drop function cannot of course capture the full complexity 
of the abundance profiles in Figure 11. 

Both models fit the CS lines despite the very different 
radial dependences of the CS abundance. The H 2 CO lines 
are better fitted by the step function models, primarily 
because they allow a high abundance over a substantial 
range of radii, where the density is high, thus providing 
stronger lines of high-excitation transitions. Similarly, the 
chemical models, even with enhanced ionization, cannot 
produce sufficient HCO + to match the observations of rare 
isotopologues. Interestingly, the chemical models do bet¬ 
ter on the HCN J = 1 —» 0 line, but neither model can 
match the peculiar J = 3 —► 2 line profile. 

In comparing models, one should bear in mind that 
the step function models were allowed 15 free parameters, 
while the chemical models enjoyed only 3, and those were 
restricted by prior knowledge. In fact, the step function 
abundances that fit best are very inconsistent with what 
we know of chemistry. The CO abundance is depleted 
inside r,„,/, while the HCO + and H 2 CO abundances in¬ 
crease; this combination is highly unlikely, especially for 
HCO + , which is a direct product of CO. The very high 
abundance of HCO + inside r, n / invoked by the step func¬ 
tion models to match the strong lines of rare isotopologues 
is very hard to produce in any chemical model. 

The evolutionary model does not include grain surface 
reactions, so it does not predict an ortho-para ratio for 
H 2 CO; we assumed a ratio of 1.7. The most likely route 
for modifications to the chemical models is to add grain 
surface reactions, but this step will effectively add many 
free parameters to the chemical model because rates for 
surface reactions are poorly known. 

On the whole, it is in fact rather remarkable that the 
line profiles from the self-consistent chemical models are 
as close as they are to the observations. With variations in 
abundances of many orders of magnitude with radius, they 
could easily have failed to match observations by orders of 
magnitude. In addition, B335 is only one source, and it 
has a rather rich spectrum for a low mass core, including 
unusually strong lines of HC ls O + . Of course, this feature 
makes B335 an attractive source to test theories, but it 
also may mean that it is atypical. 

What can explain the remaining differences between the 
models and the data? First, the Shu (1977) model of the 
infall may not be correct. There is a hint in this direction 
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in the fact that the models do not reproduce the shift of 
Vdip to the red of vlsr for optically thin lines (Fig. 10). 
This result suggests that the outer envelope is not sta¬ 
tionary, as in the Shu solution, but is also moving inward. 
Second, the outflow is not considered in these models, but 
it will affect abundances and line profiles. To include the 
outflow, we must move beyond 1-D models. Third, the 
chemical models may not yet include enough processes for 
desorption of gas from grain mantles. These could increase 
the abundances at intermediate radii. 

6.1. Off-Center Spectra 

We have focused here on spectra toward the center of 
B335. The large differences between the empirical and 
evolutionary models in the abundances of some species at 
larger radii suggest that spectra off the center of B335 may 
help in testing chemical models. Figure 12 shows a sam¬ 
ple of spectra at positions off the center predicted by the 
best-fit empirical and evolutionary models. For HCO + 
and C 18 0, we show observed spectra at these positions, 
produced by averaging spectra displaced in all directions 
with that separation in our maps. We also show predic¬ 
tions of CS spectra displaced by 60"; we no longer have 
the observed spectra at those positions, but they are pre¬ 
sented by Zhou et al. (1993). For the HCO + J = 4 —> 3 
line, the evolutionary model gives a weaker line, closer 
to the observations, than the empirical model at 15" off¬ 
set. Conversely, the evolutionary model produces stronger 
lines of C 18 0 at 30" offset than either the empirical model 
or the data. The differences are most dramatic for CS 
lines 60" away from the center. The strong decrease in CS 
abundance in the evolutionary models, seen in Fig. 11, 
produces lines that are much weaker than those of the em¬ 
pirical model. Zhou et al. (1993) detected lines stronger 
than either prediction, but the noise was fairly high. Im¬ 
proved maps of CS with good spatial resolution would be 
very helpful in further constraining models. The decrease 
in CS abundance at large radii is quite sensitive to the ex¬ 
ternal extinction (see Fig. 13 of Lee et al. 2004). Models 
with higher values of external extinction, but with all other 
parameters unchanged, do greatly increase the strength of 
CS lines at off-center positions. Maps of appropriate lines 
can help to constrain the environment of the core, while 
observations with much higher resolution can test the pre¬ 
dictions of abundance peaks at small radii (Lee, Evans, & 
Bergin 2005). 

6.2. What about the Dust? 

We noted above that the molecular line emission and 
the dust emission lead to inconsistent conclusions about 
the density distribution. To summarize, the dust emis¬ 
sion data is consistent with an inside-out collapse model 
only if the infall radius is much smaller than the molecu¬ 
lar line data indicate (Shirley et al. 2002). Attempts to 
adjust abundances to make the line profiles predicted for 
the small infall radius match the observations failed; the 
predicted line profiles are simply too narrow at the early 
times implied by the small infall radius. 

The other possibility is that the model based on the 
dust emission is wrong. While Shirley et al. performed 
an extensive set of tests, there are two possibilities that 
remain to be considered: changes in the dust opacity as a 


function of radius; and contributions to the dust emission 
from a disk. Both of these would work by adding flux at 
small radii, steepening the predicted radial intensity pro¬ 
file. For example, Young et al. (2003) found that one 
could overestimate p , the best-fit exponent in a power-law 
model by up to 0.5, if the contribution of a disk was ig¬ 
nored. This difference could change the best-fit value of 
p = 1.8 to something more compatible with the inner part 
of an inside-out collapse. 

In fact, Harvey et al. (2004) have found evidence 
for a disk in B335 and modeled multi-configuration data 
from the IRAM Plateau de Bure Interferometer with an 
envelope-disk combination. They find a good fit with a 
disk producing a flux of 21 mJy at 1.2 mm and an en¬ 
velope with a broken power law: n = 3.3 x 10 4 (r/rb)^ 1,5 
for r < rf,; and n = 3.3 x 10 4 (r/rb)~ 2 ° for r > rb . The 
best fit rb is 0.032 pc, essentially the same as our best 
fit for the infall radius, and the two broken power laws 
agree with the inside-out model for r > ri„j , and with the 
asymptotic behavior of inside-out collapse at small radii. 
However, the density just inside the infall radius has a 
slower dependence on r in the inside-out collapse solution; 
Harvey et al. (2004) note that this theoretical n(r) does 
not fit their data well. 

We have scaled up the flux of the disk in the Harvey 
et al. (2004) model (S u oc v 3 ) and modeled the 850 and 
450 /jrn data from Shirley et al. (2000) with an inside- 
out collapse model and a point source with the flux of the 
disk. The intensity profiles are not well-fitted unless the 
point source is 5-10 times stronger than in the Harvey et 
al. model. Models with disks and radial variations in dust 
opacity are needed to resolve these questions, but they are 
outside the scope of this paper. 

We conclude that the issues of disks (or more generally, 
compact structure) and radial variations in dust opacity 
introduce enough uncertainty that it is not yet possible to 
close the loop fully on the modeling of molecular lines and 
dust continuum. 

7. CONCLUSIONS 

We have assembled data on a large number of molec¬ 
ular lines toward B335 and compared those lines to pre¬ 
dictions of models. The models are of two primary types: 
those with step function chemical abundances and those 
with chemical abundances resulting from an evolutionary 
calculation (Lee et al. 2004). In both cases, the tempera¬ 
tures of dust and gas are calculated from the luminosity of 
the protostar and the density distribution under study. In 
both cases, the underlying physical model is the inside-out 
collapse model of Shu (1977). 

Both step function and evolutionary chemical abun¬ 
dances can fit most of the data, with some residual puzzles 
remaining. Both favor an infall radius around 0.03 pc, as 
was found from earlier modeling by Zhou et al. (1993) and 
Choi et al. (1995). Models with the same infall radius, but 
with densities enhanced by a factor of 5, as suggested by 
Harvey et al. (2001) do not fit as well, but they cannot 
be ruled out. Models with much smaller infall radii, as 
favored by the dust continuum modeling (Shirley et al. 
2002) do not fit the data well at all. Resolving this dis¬ 
crepancy between the conclusions of modeling dust and 
gas may require modifications to the dust models, includ- 
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ing incorporation of compact sources and changes in the 
dust opacities with radius. 

The standard chemical evolution model shows huge vari¬ 
ations in abundance as a function of radius (Fig. 11 and 
Lee et al. 2004), but still comes reasonably close to match¬ 
ing the observations. This is quite a remarkable fact. 
Changes to the initial sulfur abundance and cosmic ray 
ionization rate improve the fit to the lines, but these may 
simply be compensating for remaining unknowns in the 
chemistry. Rawlings and Yates (2001) have highlighted 
the extreme sensitivity of some abundances and line pro¬ 
files to free parameters, especially the early evolutionary 
history. Accordingly, the reader is cautioned that other 
combinations of history, dynamics, and chemical parame¬ 
ters that we have not explored might also produce reason¬ 
able fits to the data. The important point is that chemical 
models now come close enough that one can begin to look 
in detail at what might improve the match. However, this 
should be done after more than one source is compared to 
the models, as there will be variations in conditions and 
evolutionary history from source to source. 

In addition, the standard physical model of inside-out 


collapse does a remarkably good job of predicting the line 
profiles. However, there are clear hints of dynamics be¬ 
yond the standard model in the shift of the velocity of 
the self-absorption dip to higher velocities in lines requir¬ 
ing higher excitation. Models with envelopes moving in¬ 
ward may be more successful in reproducing these features. 
Spectra at positions away from the center can constrain 
other parameters, especially the surrounding radiation en¬ 
vironment. However, future work should also account for 
the non-sphericity and other effects of the outflow in this 
source. Further observations with better spatial resolution 
will be important to constrain these models, as the line 
profiles become more diagnostic of both dynamics (Choi 
2002) and chemistry (Lee et al. 2004). 

We thank the referee for perceptive comments that led 
to improvements in the paper. We thank Y. Wang, who 
collected the Haystack data and who was involved in the 
initial work on this project. We also thank Y-S. Park and 
W. Irvine for supplying data from their papers. This re¬ 
search was supported in part by NSF grants AST-9988230 
and AST-0307250 to the University of Texas at Austin. 
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Fig. 1.— Observations of the CS lines observed at IRAM (Zhou et al. 1993) (black solid histogram). The gray solid line is for the model 
with the step function abundance, and the dotted line shows the model with the chemical abundance. 
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Fig. 2. — The observations of H 2 CO are shown as black solid histograms and the models are shown as gray solid lines (step functions) and 
dotted lines (chemical abundances). The telescopes where the data were obtained are identified in each panel. The IRAM data were taken 
from Zhou et al. (1993) and the NRAO data were supplied by W. Irvine, based on data in Minh et al. (1995). The bottom two panels are 
para-H 2 13 CO. 
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Fig. 4. — The observations of HCO + are shown as black solid histograms and the models are shown as gray solid lines (step functions) and 
dotted lines (chemical abundances). The J = 1 —> 0 lines were obtained at Haystack and the other lines were obtained at the CSO. Both 
model and observations of the H 13 CO + J = 1 —> 0 line have been multiplied by 0.25 to fit them on the same scale as the 11C 13 O 1 line. 
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Fig. 6.— The observations of the CN J = 2 —*■ 1, N 2 1T 1 J = 1 —> 0 and HCN J = 1 —> 0 lines are shown as black solid histograms. For 
CN, the dotted line is a fit of Gaussians to the hyperfine components, not a true model. The fit to hyperfine components clearly does not 
reproduce the shape of the main group of lines. The dip is centered at 8.35 km s —1 , indicating that the main group is self-reversed. For N 2 H”*” 
and HCN J = 1 —> 0, the gray solid (step functions) and dotted (chemical abundances) lines are predictions of the radiative transfer model. 















T (K) Density (c 


B335: A Laboratory for Astrochemistry 



-3 - 2.5 -2 - 1.5 

log r(pc) 

The density and temperatures for gas and dust plotted as a function of radius for the standard model. 
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Fig. 9.— The observed CS lines from Zhou et al. (1993) are plotted with the line profiles predicted by various models. The step function 
models are gray and the evolutionary models are dotted. The bottom panel shows the best model with Vi n f = 0.01 pc; the second panel 
from the bottom shows the best model with ri n f = 0.03 pc (same as Fig. 1); the third panel from the bottom shows the best model with 
f'inf = 0.09 pc; and the top panel shows the best model with the ^ n / — 0.03 pc, but with densities enhanced by a factor of 5 everywhere, as 
favored by Harvey et al. (2001). 
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log r (pc) log r (pc) 



log r (pc) 


Fig. 11.— Abundance profiles. The dashed line shows the step function abundances, and the solid and dotted lines represent the abundances 
calculated from the chemical evolution model of Lee et al. (2004) in the time step for ri n f = 0.03 pc. In both chemical models, we used the 
initial sulfur abundance greater than the standard value in Lee et al. (2004) by a factor of 5. In the chemical model with dotted lines, the 
cosmic-ray ionization rate is two times greater than the standard value to give the best fit to observed line profiles (Model 6 in Table 8). The 
solid line is Model 5, which uses the standard value for ionization rate. 
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Table 1 

Observing Parameters 


Line 

V 

(GHz) 

Telescope 

T]mb 

0b 

(") 

ov 

(km s -1 ) 

Ref 

Date 

Cl 1^0 

492.1607 

CSC 

0.47 

16 

0.080 

1 

1996 Jun 

CN 2—>1 

226.874745 

cso 

0.60 

27 

0.14 

1 

1998 Jul 

C 17 0 2->l 

224.714368° 

cso 

0.81 

33 

0.17 

1 

2000 Jun 

C 18 0 2—>1 

219.560352 

cso 

0.57 

28 

0.15 

1 

1998 Jul 

C 18 0 3^2 

329.3305453 

cso 

0.82 

26 

0.10 

1 

2000 Jul 

HCO+ l->0 

89.188512 

Haystack 

0.12 

25 

0.10 

1 

1994 Jun 

H 13 CO+ l->0 

86.754330 

Haystack 

0.12 

25 

0.10 

1 

1994 Jun 

HCO+ 3-»2 

267.557619 

CSO 

0.65 

26 

0.18 

1 

1995 Mar 

H 13 CO+ 3—>2 

260.255339 

CSO 

0.65 

26 

0.18 

1 

1995 Mar 

HC 18 0+ 3—>2 

255.47940 

CSO 

0.65 

26 

0.18 

1 

1995 Mar 

HCO+ 4-+3 

356.734288 

CSO 

0.61 

20 

0.14 

1 

1995 Mar 

DCO+ 3->2 

216.112604 

CSO 

0.57 

28 

0.16 

1 

1998 Jul 

HCN l->0 

89.635847 

TRAO 

0.40 

61 

0.068 

2 

1997 

HCN 3—>2 

265.8864343 

CSO 

0.65 

23 

0.15 

1 

1996 Jun 

H 13 CN 3^2 

259.011814 

CSO 

0.65 

23 

0.15 

1 

1996 Jun 

HNC 3—»2 

271.981142 

CSO 

0.62 

22 

0.11 

1 

1996 Jun 

HNC 4—>3 

362.630303 

CSO 

0.53 

19 

1.62 

1 

1997 Jul 

HC 3 N 24^23 

218.324788 

CSO 

0.56 

28 

0.18 

1 

1996 Jun 

N 2 H+ l->0 

93.176258 b 

Haystack 

0.12 

25 

0.10 

1 

1994 Jun 

N 2 H+ 3^2 

279.511757° 

CSO 

0.56 

22 

0.12 

1 

1996 Oct 

N 2 D+ 3—>2 

231.321775 

CSO 

0.73 

32 

0.21 

1 

2001 Jul 

para-H 2 13 CO loi - 0oo 

71.02478 

NRAO 

0.95 

89 

0.206 

3 


H 2 13 CO 2i 2 — In 

137.44996 

NRAO 

0.72 

42 

0.1065 

3 


H 2 CO 2i 2 — In 

140.839518 

IRAM 

0.68 

17 

0.083 

4 


para-H 2 13 CO 2 q 2 — lgi 

141.98375 

NRAO 

0.72 

42 

0.103 

3 


H 2 13 CO 2n — lio 

146.63569 

NRAO 

0.72 

42 

0.0998 

3 


H 2 CO 3i 2 — 2n 

225.697772 

IRAM 

0.50 

12 

0.066 

4 


H 2 CO 3i 2 — 2n 

225.697787 

CSO 

0.65 

27 

0.127 

1 

1996 Jun 

para-H 2 CO 3 03 - 2 02 

218.222186 

CSO 

0.65 

28 

0.131 

1 

1996 Jun 

H 2 CO 5is - 4 14 

351.768645 

CSO 

0.53 

20 

0.083 

1 

1997 Jun 

para-H 2 CO 5 0 5 - 4 04 

362.3530480 

CSO 

0.53 

19 

0.085 

1 

1997 Jun 

para-H 2 CO 5 23 - 4 22 

365.3634280 

CSO 

0.53 

19 

0.085 

1 

1997 Jun 


Note. — (a) Reference frequency for the hyperfine shifts in Ladd et al. (1998); (b) For the isolated hyperfine component (Lee et al. 2001); 
(c) Reference frequency for the hyperfine shifts in Caselli et al. (2002); (1) This paper; (2) Park et al. 1999; (3) Minh et al. 1995; (4) Zhou et 
al. 1993. 



B335: A Laboratory for Astrochemistry 


23 


Table 2 

Observational Results 


Molecule 

Line 



J T\dv 
(K km s -1 ) 

T* a 

(K) 

vlsr 
(km s -1 ) 

Av 
(km s -1 ) 

Cl 

.J = \ -, 

0 


0.80(0.04) 

1.51(0.18) 

8.29(0.01) 

0.50(0.03) 

CN fc,c 

J = 2 -> 

1 


0.90(0.02) 

0.43(0.02) 

8.35(0.04) 

0.62(0.07) 

c 17 o° 

J = 2 —» 

1 


0.46(0.07) 

0.57(0.07) 

8.39(0.02) 

0.49(0.07) 

c 18 o 

J = 2-> 

1 


0.63(0.42) 

1.10(0.10) 

8.27(0.02) 

0.57(0.04) 

c 18 o 

J = 3-> 

2 


0.78(0.09) 

2.15(0.26) 

8.33(0.04) 

0.64(0.08) 

HCO +h 

J = 1 -> 

0 


0.71(0.05) 

1.15(0.11) 

8.35(0.03) 

0.61(0.10) 

h 13 co+ 

J = 1 -+ 

0 


0.24(0.06) 

0.33(0.02) 

8.23(0.04) 

0.58(0.03) 

HCO+ 6 

J = 3 -> 

2 


3.13(0.05) 

3.32(0.08) 

8.46(0.04) 

0.94(0.16) 

h 13 co+ 

J = 3 —> 

2 


0.46(0.01) 

0.76(0.03) 

8.25(0.01) 

0.57(0.02) 

hc 18 o+ 

J = 3 —> 

2 


0.05(0.01) 

0.09(0.01) 

8.26(0.04) 

0.58(0.10) 

HCO+ 6 

J = 4-» 

3 


2.27(0.15) 

2.35(0.06) 

8.55(0.03) 

0.97(0.12) 

DCO+ 

J = 3-, 

2 


0.78(0.03) 

1.32(0.04) 

8.35(0.01) 

0.55(0.02) 

HCN fc,c 

J = 1 -> 

0 


0.59(0.10) 

0.35(0.05) 

8.39(0.01) 

0.70(0.10) 

HCN h 

J = 3^> 

2 


0.81(0.03) 

0.80(0.02) 

8.40(0.20) 

1.01(0.04) 

h 13 cn 

J = 3 -» 

2 



< 0.1 



HNC 

J = 3-> 

2 


0.22(0.01) 

0.49(0.03) 

8.16(0.01) 

0.42(0.03) 

HNC 

J = 4-> 

3 


0.39(0.05) 

0.16(0.03) 

8.33(0.16) 

2.3(0.4) 

hc 3 n 

J = 24 - 

v 23 



< 0.08 



N 2 H +d 

J = 1 

0 


0.13(0.01) 

0.25(0.02) 

8.33(0.01) 

0.47(0.02) 

N 2 H+ b ' c 

J = 3-> 

2 


0.95(0.08) 

1.00(0.04) 

8.38(0.03) 

0.38(0.04) 

N 2 D+ a 

J = 3 — > 

2 


0.32(0.04) 

0.38(0.04) 

8.36(0.02) 

0.31(0.05) 

h 2 co 

Jk- x k +1 

= 3i 2 - 

-> 2n 

0.55(0.01) 

0.63(0.02) 

8.26(0.01) 

0.82(0.02) 

para-H 2 CO 

Jk^k +1 

= 3o3 - 

-+ 2 02 

0.45(0.01) 

0.59(0.03) 

8.28(0.01) 

0.71(0.03) 

H 2 CO 

Jk-!K +1 

= 5is - 

4i4 

0.49(0.03) 

0.57(0.07) 

8.30(0.02) 

0.80(0.07) 

para-H 2 CO 

Jk^ x k +1 

= 5o5 - 

-+ 4 0 4 

0.40(0.06) 

0.31(0.11) 

8.53(0.08) 

1.22(0.23) 

para-H 2 CO 

Jk_!K +1 

= 5 23 - 

-+ 4 22 


< 0.2 




Note. — (a) Hyperfine structure: f T^dv refers to area under all lines; TJ refers to main peak of blended components; vlsr and Ac are 
from fit to blended components; (b) Double-peaked: f T^dv is for total area under both peaks; TjJ refers to strongest peak; vlsr refers to 
dip; Av is integrated intensity divided by peak TJ. (c) Hyperfine structure: f T\dv refers to area under all lines; At) determined from isolated 
hyperfine component, (d) All entries refer to the isolated hyperfine component. 


Table 3 

Standard Physical Model 


Type a 

(km s -1 ) 

Tinf 

(pc) 

r out 

(pc) 

L 

(Le) 

Ay 

(mag) 

Go C 

(s 3 ) 

b 

(km s -1 ) 

n e 

(cm -3 ) 

X(CO) 

Shu 0.23 

0.03 

0.15 

4.5 

1.3 

“or 3 x io -iY 

0.12 

1 x 10 -3 

7.4 x 10 -5 


oT^o 

70 540 


i 8 0/ iV 0 

3.5 


Table 4 
Isotope Ratios 
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Table 5 

Aggregated Hyperfine Shifts 3, and Strengths 15 


C iY 0 

2^1 

HCN 

l->0 

HCN 

3^2 

n 2 h+ 

1-V) 

n 2 h+ 

3^2 

Av 

Ti 

Av 

n 

Av 

n 

Av 

Ti 

Av 

n 


1.157 

0.040 

4.849 

0.333 

1.749 

0.037 

6.936 

0.037 

2.015 

0.017 

0.431 

0.122 

0.000 

0.556 

0.303 

0.200 

5.984 

0.185 

0.669 

0.015 

0.241 

0.571 

-7.072 

0.111 

-0.030 

0.725 

5.545 

0.111 

0.416 

0.084 

-0.526 

0.093 



-0.611 

0.001 

0.956 

0.185 

0.266 

0.094 

-0.926 

0.016 



-2.348 

0.037 

0.000 

0.259 

0.076 

0.089 

-1.073 

0.095 





-0.611 

0.111 

-0.073 

0.615 

-1.203 

0.062 





-8.006 

0.111 

-0.601 

0.010 









-2.644 

0.011 









-2.773 

0.014 


Note. — (a) Shift in km s 1 relative to the assumed central velocity of 8.30 km s 1 ; (b) Relative strength (r I: ) normalized so that 


Table 6 

Step Function Abundances 


Species 

X{r > r inf ) 

X(r < ri n f) 

AD(Step) 

AD(Chem) 

CS 

6.0 x 10~ 9 

6.0 x 10 -9 

6.88 

6.78 

C 18 0 

7.4 x 10~ 8 

2.5 x 10" 8 

3.00 

3.22 

HCO+ 

3.5 x 10“ 9 

3.5 x 10" 8 

3.95 

4.02 

DCO+ 

6.0 x 10" 11 

6.0 x 10“ 10 

1.87 

3.27 

n 2 h+ 

3.0 x 10~ 10 

3.0 x 10" 9 

5.15 

6.14 

HCN 

1.0 x 10~ 8 

i.o x nr 9 

5.43 

4.73 

H 2 CO 

7.0 x 10" 10 

7.0 x 10" 9 

2.64 

3.22 

para-H 2 CO 

4.0 x 10" 10 

4.0 x 10" 9 




Table 7 

Variations in the Physical Model 


r inf 

Density factor 

(AD) a Chem Mod. 

(. AD) a 

0.01 

1.0 

4.62 

6 

4.93 

0.03 

1.0 

3.75 

6 

4.15 

0.09 

1.0 

5.38 

6 

4.95 

0.03 

5.0 

4.02 



The mean value of AD for all 

species. 
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Table 8 

Chemical Models 


Model 

Type 

Dust 

Surface 

X(S) 

c 

(s- 1 ) 

{AD) a 

Model 1 

Step 




3.75 

Model 2 

Chem 

Si0 2 

4 x 10" 8 

3.0 x 10~ 17 

5.05 

Model 3 

Chem 

CO 

4 x 10" 8 

3.0 x 10- 17 

4.84 

Model 4 

Chem 

H 2 0 

4 x 10" 8 

3.0 x 10- 17 

5.89 

Model 5 

Chem 

Si0 2 

2 x 10" 7 

3.0 x 10- 17 

4.25 

Model 6 

Chem 

Si0 2 

2 x 10" 7 

6.0 x 10- 17 

4.15 

Model 7 

Chem 

Si0 2 

2 x 10" 7 

1.5 x 10~ 16 

4.23 


“The mean value of AD for all species. 



